d-wave Superfluid with Gapless Edges in a Cold Atom Trap 
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We consider a strongly repulsive fermionic gas in a two-dimensional optical lattice confined by 
a harmonic trapping potential. To address the strongly repulsive regime, we consider the t — J 
Hamiltonian. The presence of the harmonic trapping potential enables the stabilization of coexisting 
and competing phases. In particular, at low temperatures, this allows the realization of a d-wave 
superfluid region surrounded by purely (gapless) normal edges. Solving the Bogoliubov-de Gennes 
equations and comparing with the local density approximation, we show that the proximity to the 
Mott insulator is revealed by a downturn of the Fermi liquid order parameter at the center of the 
trap where the d-wave gap has a maximum. The density profile evolves linearly with distance. 



Ultracold atoms in optical lattices are ideal quan- 
tum simulators of complex many-body Hamiltonians that 
arise in condensed-matter systems [TH3]- They embody 
very clean systems which can be tuned in a very precise 
and controlled manner from the weak to the strong cou- 
pling limit. This enables one to investigate a plethora 
of interesting phenomena ranging from the dynamics 
of strongly correlated bosons [H [5] and fermions [S [7] 
to quantum magnetism [8J. The pace of experimental 
progress is quite impressive. Bunching and antibunch- 
ing effects in the density-density correlations were found 
for bosons and fermions [5]. In particular, the fermionic 
Hubbard model has been realized for repulsive and at- 
tractive interactions. Fingerprints of the Mott state 
[TU1 HTj and s-wave superfluidity [H] have already been 
observed in experiments. The challenge remains to ac- 
cess the Neel phase [13] in order to reveal the presence of 
high-temperature d-wave superfluidity close to half-filling 
|14j . Spin fluctuations, which are predicted to be the glue 
for d-wave superfluidity in the repulsive Hubbard model 
|15H18j . have been studied experimentally in the context 
of BEC-BCS crossover fTS]. In this Letter, we theoreti- 
cally address the effect of the harmonic potential, which 
originates from the Gaussian profile of the laser beams 
generating the trap, on the d-wave superfluid phase of the 
Hubbard model. We explore the strongly repulsive limit 
which is realized in high-T c superconductors |15H18j . 

We consider relatively small fillings such that only a 
superfluid cloud is left in the center of the trap and anti- 
ferromagnetism is hindered by the motion of atoms |20j . 
We investigate the coexistence between d-wave superflu- 
idity and a normal phase at the boundaries applying an 
effective theory of a doped Mott insulator. More pre- 
cisely, we start from the t — J Hamiltonian and apply a 
low-energy (superfluid) theory that allows us to describe 
the proximity to the Mott insulator [2"TH2"3"] . 
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where c ia creates a particle with spin a at the i'th lat- 
tice site, t, the tunneling amplitude, is chosen to be spin- 
and direction-independent, Vi denotes an isotropic and 
spin-independent external trapping potential, V(x,y) = 
(m/2)uj 2 (x 2 + y 2 ) evaluated at the i'th lattice site, 
and (ij) denotes nearest neighbor pairs. The last 
term describes the super-exchange process where Si = 

^/2J2cra' HtT a ao'Cia' ls tne spin-operator and J = At 2 /U. 
The number of fermions per site is assumed to be less 
than one. The ground state properties can be ana- 
lyzed using a projected Bardeen-Cooper-Schrieffer (BCS) 
state taking into account that double occupancy is for- 
bidden for large interactions, Pq |BCS) where Pq — 

n, (i - n if %) us mi mi- 

We account for the strong repulsion by projecting onto 
the states with site-population less than or equal to 
unity which is done by introducing the site-dependent 
Gutzwiller factors [21] : 

fl « = v& = ^7(1 + ^)^/(1 + ^) (2) 

= JgJi = V4/(1 + ^)V4/(1 + ^) 2 . 

where we denote the deviation from half-filling 6 l — 
1 — (hi) with (hi) — J2 a (fticr) representing the total 
number of particles on site i. This is justified by choos- 
ing a small curvature of the external trapping potential 
compared to the hopping amplitude. This procedure is 
shown to be in good agreement with variational Monte 
Carlo calculations for the projected d-wave state [25] . 
The renormalized Hamiltonian then takes the form 

H = -t^2 \9l 3 c\aCjo + h.c. 

+ J2v t clc ta + jJ2gl 3 s i -s 1 , (3) 
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where gl J works as to enhance spin-spin correlations since 
gl > 1 and g\ 3 suppresses the tunneling term since g\ 3 < 
1 and goes to zero when the doping goes to zero. This 
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resembles the suppression of the kinetic energy in the 
Mott insulating state. 

In order to solve this many-body Hamiltonian we in- 
troduce the Fermi liquid order parameter, Xiji an( l the 
pairing order parameter, A^-, with the following averages 

EH 



Xij 
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where = 1 = —e^ and e-f-j- = = e^. At a gen- 
eral level, spin fluctuations not too far from half-filling 
will favor d-wave superconductivity. We then obtain the 
mean-field Hamiltonian 
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where we now work in the grand-canonical ensemble and 
therefore introduce the chemical potential \i since the 
number of particles is no longer conserved. We have in 
the above expression assumed that the potential is suf- 
ficiently slowly varying so that (cj-fCjj,) = (cj-j-Cjj,). Fur- 
thermore, we have exploited the fact that the additional 
term to the tunneling amplitude is spin-independent 
(cLcj-f) = (c\^Cji) under our present assumptions. 

Below, we show that the pairing order parameter de- 
pends strongly on the position in the trap. It converges 
to the untrapped case in the centre of the trap where it 
is maximal, which reflects the emergence of superfluidity 
induced by the pairing term. We confirm the co-existence 
of a d-wave superfluid domain and a gapless edge region. 

Since the trap breaks the discrete translational invari- 
ance, we first find the gap by solving the Bogoliubov- 
de Gennes (BdG) equations. We diagonalize the renor- 
malized mean-field Hamiltonian by making a unitary 
Bogoliubov-Valatin transformation of the creation- and 
annihilation operators, expanding them on a complete 
basis of quasi-particle modes annihilated by the fermionic 
mode operators j V(r [25] 
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Here the sums are performed solely over positive energies 
and it^! is the amplitude of destroying a quasi-particle 
with spin up at the i'th lattice site. The self-consistent 
equations of the mean-field parameters take the following 
form in terms of the quasi-particle amplitudes 



xn = ^ E [<Xt/ ( E m) + <4</ (-Ent) 
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where f(E vt ) = ftJtV) = l/(™p(E vt /T) + 1) is the 
thermal occupation of the quantum state with energy 
E^f. In the above expression the duality of the spin- up 
and spin-down solutions has been exploited to transform 
the sum into a sum over positive as well as negative en- 
ergies for spin-up particles only. 

Assuming a small curvature of the harmonic trapping 
potential, it is possible to approximate the potential with 
a constant in the vicinity of each lattice site and replace 
the eigenstates with plane wave states. Assuming a large 
number of particles in each cell of constant potential, we 
define a local Fermi sea. This consists in defining a local 
chemical potential in each cell fj, (r) = fj, — V(x,y) and 
letting 
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where u k {r) and v k (r) are solutions of the homogenous 
system with /z(r) = [i— V(x, y) [2Z] and Nl is the number 
of lattice sites. Below, we introduce the local variables 
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where a denotes the nearest neighbors, Njy the number of 
nearest neighbors, and adx — oi-dx = —ot c i y = —ot-d y 
1 

Within this local density approximation (LDA) the 
self-consistent equations take the form 
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FIG. 1: (Color online) A diagonal cross-section of A(r) of the 
trapped system with 76 x 76 lattice sites, as a function of the 
distance from the center of the trap ro. t/J — 5, fi/J = —0.1, 
and ui = yf ma 2 / Jlo — 0.16. BdG (circles) and LDA (crosses). 

with 



Ek(x,y) =\/£,k{x,y) 2 + A k (x,y) 2 
£.k{x, y) = - (2tg t (x, y) + x(x, y)) [cos(fc x a) + cos(k y a)} 

Afc (x, y) =A(x, y) [cos(k x a) - cos(k y a)] . (13) 

We have neglected the variation in chemical potentials 
on neighboring lattice sites and approximated the renor- 
malization factors with their value on the i'th lattice site, 
for example g s (x, y) = 4/(1 + 6(x, y)) 2 , which is justified 
when the curvature of the harmonic potential is small 
compared to the tunneling amplitude. In the absence of 
the harmonic trap, these equations and their solutions 
agree with well-known results; see for example Ref. |23j . 

We now compare the local density approximation cal- 
culation to the BdG calculation of the site-dependent val- 
ues of Aj, Xi an d (^i)i n °t e that the tilde values mean 
that the mean-field order parameters are expressed in 
units of 3/4<?" J. For large interactions, we check that the 
order parameters are controlled by the Anderson super- 
exchange J. We expect that the critical distance at which 
A, will be zero is given by the position with a doping cor- 
responding to the critical "doping" from the homogenous 
calculation which lies around S = 0.35 [23 if one assumes 
the LDA is correct. Remarkably, we obtain a very good 
quantitative agreement between the LDA approximation 
and the BdG calculation as shown in Figs. 1 and 2. We 
corroborate the coexistence between <i-wave superfluid- 
ity and normal normal edges in accordance with results 
from a calculation on a similar system sufficiently away 
from half- filling in the weakly interacting limit [20] ■ In 
addition, we find that in the center of the trap, since the 
fermion density reaches one, there is a reminiscence of 
the Mott insulating state which results in a downturn of 
the Fermi liquid order parameter \i- 

In the following we consider the variation of the dop- 
ing close to the boundaries of the density profile de- 
fined as the circle in the ccy-plane with radius R where 
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FIG. 2: (Color online) A diagonal cross-section of x(r) of the 
trapped system with 76 x 76 lattice sites, as a function of the 
distance from the center of the trap ro. t/J = 5, fi/J = — 0.1, 
and ui = y/ma 2 / Jlu — 0.16. BdG (circles) and LDA (crosses). 

5(x,y) reaches unity in the local density approximation. 
This circle is only well-defined in the continuous sys- 
tem but becomes sufficiently well-defined for large lat- 
tices. In the vicinity of the density boundaries the pair- 
ing order parameter has vanished, x( x ,y) is so small 
that we ignore it and at the boundaries the chemical po- 
tential n(R) is equal to —At. This means that all the 
terms in the expression for 5(R) contribute with a posi- 
tive sign since £k{R) = —2t (cos(k x a) + cos(k y a)) + At > 

for all quasi-momenta and S(R) is therefore maxi- 
mal and unity. In a small region close to 5(R) = 

1 the terms in the sum which contribute to a low- 
ering of the doping concentration will form a circle 
in fc-space. Hence, we make a Taylor expansion of 
the cosine terms in fc-space at site i with coordinates 
{x, y) to find the radius ka corresponding to the change 
of sign in Ck(x,y)- -2tg t (x, y)(cos(k x a) + cos(k y a)) - 
H(x,y) w -2tg t (x,y)(2 - 1/2 [(k x a) 2 + (k y a) 2 ] ) - 
fi(x,y) = -2tg t (x,y)(2 - \/2{ka) 2 ) - ^{x,y) = which 



corresponds to a radius of (ka) 2 
ing will then be given by 
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with [i(x,y) = ji — V(x, y) the local chemical potential 
when realizing that one must subtract the area corre- 
sponding to the negative contribution twice. 
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FIG. 3: (Color online) A diagonal cross-section of (n(r)} of the 
trapped system with 76 x 76 lattice sites, as a function of the 
distance from the center of the trap ro. t/J — 5, fi/ J = —0.1, 
and lu — \J ma 2 / Juj = 0.16. BdG (circles), LDA (crosses) 
and analytical (dashed). 

One can use the expression for the chemical potential 
at the boundaries n{R) = fi— l/2mw 2 fl 2 = —At to obtain 
the associated radius R = y / 2(4t + /i) / (mw 2 ), which 
corresponds to a radius of R = 39.5a when t/J = 5, 
T/J = 0, n/J = —0.1, and Q = 0.16 which is shown to 
be correct in Fig. [3] 

Due to the isotropy of the harmonic potential the ex- 
pression for the doping is assumed to have no angular- 
dependence as in the continuous case and the expression 
for the density can be Taylor expanded around R to first 



order in x to give n(x,0) = — 0.024(x/a — 39.5). 

A qualitative explanation of A, is given at finite tem- 
perature. We observe that the height and the width of Aj 
decreases with increasing temperature whereas the shape 
does not change significantly. 

To summarize, we have shown that a harmonic trap- 
ping potential can simulate the effect of inhomogeneous 
doping allowing to stabilize a novel d-wave superfluid 
phase for fermions with gapless edge states. In the con- 
text of strongly repulsive fermions, we have found that 
at the center of the trap there is a reminiscence of the 
Mott insulating state which is revealed by a downturn of 
the Fermi liquid order parameter whereas the d-wave gap 
is maximal. The proximity to the Mott insulating state 
at the center of the trap should affect the double occu- 
pancy. The d-wave symmetry of the pairing can be de- 
tected through the bunching which is maximal along the 
x-and y-directions and minimal along x = ±y. The pair- 
ing should be maximal around the Fermi surface. Phase 
sensitive measurements could also be performed [28] . We 
found a linear profile of the fermion density close to the 
boundaries. When decreasing the density, one might ob- 
serve magnetic real-space shell structures competing with 
the superfluid phase [2U]. For dipolar fermions, d-wave 
bond order solidity might also occur at half-filling |29j . 
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